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Abstract 

The variational grid generation method is a powerful tool for generat- 
ing structured convex grids on irregular simply connected domains whose 
boundary is a polygonal Jordan curve. Several examples that show the 
accuracy of a difference approximation to the solution of a Poisson equa- 
tion using these kind of structured grids have been recently reported. In 
this paper, we compare the accuracy of the numerical solution calculated 
by applying those structured grids with finite differences against the the 
solution obtained with Delaunay-like triangulations on irregular regions. 

1 Introduction 

For the numerical solution of the Poisson equation on irregular domains, the 
use of finite differences and finite elements with triangulations obtained by sub- 
dividing each grid cell of a structured grid along a diagonal has been addressed 
in Ref. [8^, proving that the finite difference approach on such grids is accu- 
rate enough. However, since structured grids often have some elongated cells, a 
natural question that arises in this context is whether this numerical solution is 
more accurate than the solution obtained by using finite elements on a standard 
Delaunay-like triangulation. 

In this paper, we compare the accuracy of the numerical solution for this problem 
using finite differences in structured grids, and finite elements on Delaunay-like 
triangulations. We are specifically interested in irregular boundaries, since their 
geometry is closer to a realistic domain, for instance, a lake. 
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In order to generate the structured convex grids, a variational method was used. 
The latter consists of minimizing an appropriate functional [6j IH [TOl [13] . Area 
and harmonic functionals can be used for gridding a wide variety of simple con- 
nected domains in the plane [H [SI HI [TTJ [TH [T7] , whose boundaries are closed 
polygonal Jordan curves with positive orientation. 

If m and n represent the "vertical" and "horizontal" numbers of points of the 
"sides" , then the boundary is the positively oriented polygonal curve 7 of ver- 
tices V = {vi, ■ ■ ■ , V2[m+n~2)} 1 ^^^d it defines the typical domain n. 
A doubly indexed set 

G = {Pij\l <i < m,l<j<n} 

of points of the plane with the fixed boundary positions given by y is a logically 
rectangular structured grid with quadrilateral elements for f2, of order m x n. 
A grid G is convex if and only if each one of the (to — l)(n — 1) quadrilaterals 
(or cells) Cij of vertices {Pij, Pi+ij , Pij+i, Pi+ij+i}, 1 < z < m, 1 < j < n, is 
convex and non-degenerate (See fig. [1]). 

The basis for the direct optimization method, as developed by Ivanenko et al. 




Figure 1: Structured grid generated by UNAMALLA 
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[S], is the minimization a suitable function of the form 



m— 1 n— 1 



(1) 



where Cij is the {i^jy^ grid ceh and / is a function of its vertices; the problem 
is to find the coordinates of the interior points of the grid G. The functional 
used to generate the structured grids of the numerical tests, as implemented in 
UNAMALL A [18] , was the adaptive linear convex combination of the area func- 
tional Sui{G) described in Ref. [5], and the length functional L{G) with weight 



The parameter w, which is a scale factor, can be updated in such a way that 
in a finite number of updates the combined functional attain its minima within 
the set of convex grids for il if the latter is nonempty. Further properties of 
this functional, as well as the algorithm for updating its parameter has been 
reported in Ref. [5] and Ref. |S]. 

The triangulations we considered for this paper are those generated by DistMesh 
[T^ , which is based on the physical analogy between a simplex mesh and a truss 
structure, where the meshpoints are nodes of the truss. It generates an initial 
Delaunay triangulation, then assumes an appropriate force-displacement func- 
tion for the bars in the truss at each iteration, and finally solves for equilibrium 
(See fig. In order to have grids with a similar number of elements, we used 
two initializations: a) the initial edge length for DistMesh was set in propor- 
tion to half the average diagonal length of the structured grids, b) a variation 
of DistMesh was designed for which the initial points inside the region are the 
inner nodes of the corresponding structured grid. 

It must be noted that for very irregular boundaries, DistMesh might produce a 
few triangular elements along the boundary which do not satisfy the Delaunay 
condition. 



2 Finite difference approximation 

Standard difference schemes can be generalized by considering a finite set of 
nodes pi,p2, ...,pk, for which it is required to find coefficients Fi, F/j such 

that 



As it is well known, the F values can be calculated with ease in regular regions. 
However, despite the basic idea is quite simple, the application to Taylor's The- 
orem leads to more complicated schemes on irregular regions; up to our knowl- 
edge, there are few efficient schemes for such kind of regions. 
The calculation of these coefficients has been studied by Tinoco ei aZ [Ij, and 
Shashkov [M]. We make use of the second order scheme developed for the 



(7 = 0.5 (See Ref. 0). 
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Figure 2: Triangulation generated by DistMesh 

boundary value problem 

-ViK{x,y)Vu{x,y)) = f{x,y) (3) 

u(a;,y)|an = g{x,y), 

with non-diagonal matrices K(x,y) studied in Ref. |14j . which is based on the 
method of support-operators and has the advantage of providing explicit for- 
mulas for the r coefficientf0 . 

In the numerical experiments, we selected the 3x3 subgrid defined for the 
nodes x{i — 1 : i + I, j — I : j + 1), y{i — 1 : i + 1, j — 1 : j + 1) to approximate 
—'S/{K{x, y)Vu{x, y)) = f{x, y) at the inner grid node {x{i, j),y{i, j)). As in the 
rectangular case, an algebraic system of equations is obtained from discretiza- 
tion, which becomes sparse as m and n increase. 



^This scheme is second order according to the grid norm defined by equation Q. 
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3 Finite elements approximation 

Let Ne be the number of triangular elements in a grid. Galerkin's approach to 
the solution of equation ([3]) is given by the combination 

Ne 

u{x,y) Ki^U^(f)i (4) 

1=1 

of trial-test pyramid functions whose faces are defined on a triangle by 

(l){x,y) ^ A + Bx + Cy. (5) 
This selection yields the weak formulation 

Ne . „ 

y^U, <Vc^j,K{x,y)V<i>,>dA= <P,fdA, j-l,---,iVe, (6) 
j^]^ JO. Jn 

where < •, • > is the canonical inner product (See Ref. |15)V 
There are several other possible choices of elements. However, as it is well 
known, the use of pyramids on triangulations allows a very simple second order 
approximation to the solution of equation ([3]). 



4 Numerical tests 

For the numerical tests, we have selected 9 polygonal regions, most of them 
approximations to real geographical locations: they will be denoted as Dome 
(dom). Great Britain (eng), Havana bay (hab), M19 (ml9), Mexico (mex). Plow 
(plo). Swan (swa), Ucha (uch) and Michoacan (mic). They are shown in figure 
3. 

Scaling these boundaries in order to lie in [0, 1] x [0, 1], convex grids with 21, 
41 and 81 points per side were generated by minimizing the 1/2(S'^(G') + L{G)) 
functional in UNAMALLA with default parameters. The resulting structured 
grids were used with Shashkov's finite difference schemes [13]. As mentioned 
before, they were also used as initial data for some triangulations. 
DistMesh was used to triangulate the same test polygonal boundaries with de- 
fault parameters, setting the bar length as half the average diagonal length 
calculated in the corresponding structured grids (denoted as DistMesha). A 
variation of DistMesh was also designed, for which the initial points inside the 
region are the inner nodes of the corresponding structured grid: these grids are 
denoted as DistMeshb grids. 

The number of elements and inner nodes in each grid can be seen in in table 
[TJ The column N gives twice the number of grid cellfl and the columns Na 

■^This number is considered because the number of triangles in a triangulation obtained by 
subdividing each grid cell along a diagonal of a structured grid is twice the number of cells. 
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Dome Great Britain Havana bay 




M19 Mexico Michoacan 




Plow Swan Ucha 

Figure 3: Test regions. 



and Nh the number of triangular elements in the DistMesha and DistAIeshb 
triangulations, respectively. The columns labeled Nu, Nua and Nuh represent 
the corresponding number of inner nodes, which is equal to the number of un- 
knowns in the approximation. 

The resulting algebraic systems for finite differences and elements are sparse 
due to the discretization. However, the systems for finite differences are block- 
tridiagonal, so they can be solved with a number of algorithms in a very efficient 
way; this is a clear advantage of the double index in a structured grid. For the 
tests, Gauss-Seidel Method was used. 

On the other hand, the matrices in the finite element systems have no specific 
structure and require a more complicated data structure in order to solve them 
efficiently. For the tests, they were solved using a sparse Gaussian elimination 
routine. 

The following values for u and K were selected (See Ref. [HI): 
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Table 1: Number of grid elements and inner nodes 



Region 


oize 


iVi 


iV2 




AT 11 
iV U 


AT 11 




dom 


21 


800 


1540 


795 


361 


693 


356 




41 


3200 


6064 


3189 


1521 


2884 


1509 




81 


12800 


24150 


12783 


6241 


11785 


6218 


eng 


21 


800 


841 


771 


361 


355 


337 




41 


3200 


3463 


3111 


1521 


1581 


1440 




81 


12800 


13918 


12614 


6241 


6622 


6074 


hab 


21 


800 


1261 


767 


361 


559 


333 




41 


3200 


5048 


3115 


1521 


2370 


1446 




81 


12800 


20313 


12605 


6241 


9845 


6070 


ml9 


21 


800 


1420 


785 


361 


626 


350 




41 


3200 


5717 


3157 


1521 


2696 


1483 




81 


12800 


22772 


12705 


6241 


11067 


6156 


mex 


21 


800 


477 


686 


361 


183 


273 




41 


3200 


1781 


2911 


1521 


771 


1288 




81 


12800 


7036 


11918 


6241 


3301 


5519 


mic 


21 


800 


1432 


796 


361 


645 


356 




41 


3200 


5760 


3183 


1521 


2747 


1506 




81 


12800 


22812 


12694 


6241 


11096 


6136 


plo 


21 


800 


863 


748 


361 


363 


310 




41 


3200 


3335 


3076 


1521 


1535 


1399 




81 


12800 


13115 


12537 


6241 


6303 


5983 


swa 


21 


800 


1410 


796 


361 


634 


357 




41 


3200 


5495 


3187 


1521 


2612 


1508 




81 


12800 


21737 


12766 


6241 


10580 


6205 


uch 


21 


800 


1088 


794 


361 


483 


355 




41 


3200 


4252 


3166 


1521 


1976 


1489 




81 


12800 


17160 


12729 


6241 


8288 


6172 
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1. First problem. 



Kix,y)^[ Q ^ ), u = 2exp(2x + 2/). 



2. Second problem, 
with 

and 



K{x,y) = P^DP, 



p ^ ( cos(|) sm(f ) 



-sm(f) cos(f) 

1 + 2x2+2/2 Q 
1 + + 2?/2 

u ~ sin(7rx) sin(7ry). 



3. Third problem, 
with 

and 

D = 



K{x,y) = P^DP, 

cos(f) sinil) 
-sm(f) cos(f) 

1 + 2x2 + y2 _|_ y5 

1 + a;2 + 2?/2 

u — sin(7ra::) sin(7r?/). 

Function / was chosen in such a way that u was the exact solution in every 
case. 

The II ■ II 2 error norm for the tests is summarized in tables [U [3] and |H it was 
calculated as a grid function 



\\u^Uh = JJ2iu.-U0'A. , (7) 



where u and U are the approximated and the exact solution calculated at the 
i'^-element, and Ai is the area of the element. The approximation correspond- 
ing to label Structured are those of the second order finite differences men- 
tioned in section [21 the approximations for the columns labeled DistMesha, 
and DistMeshb were calculated with finite elements using the pyramid trial- 
test approximation described in section [3l The empirical orders O, Oa and Ob 
between two consecutive grid orders were calculated according to the formula 
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Table 2: Quadratic error for problem 1 



Region 


oize 


S tVllCtllTCd 






UlSblVi CSfla 






UlSblVi CSfl}) 




dom 


21 


4.59E-03 






5.77E-03 






6.25E-03 






41 


1.22E-03 


1 


.98 


8.79E-04 


2, 


.81 


1.47E-03 


2.16 




81 


1.46E-04 


3 


.12 


1.51E-04 


2 


.58 


2.70E-04 


2.49 


eng 


21 


2.58E-03 






1.72E-03 






2.44E-03 






41 


4.27E-04 


2, 


.69 


6.95E-05 


4, 


.79 


4.75E-04 


2.45 




81 


1.17E-04 


1 


.90 


1.03E-04 


-0, 


.57 


2.86E-04 


0.75 


hab 


21 


7.53E-03 






1.26E-03 






1.98E-03 






41 


1.53E-03 


2, 


.38 


1.69E-04 


3 


.00 


5.85E-04 


1.82 




81 


3.97E-04 


1 


.99 


6.11E-05 


1 


.50 


1.77E-04 


1.76 


mlQ 


21 


5.90E-03 






2.02E-03 






1.94E-02 






41 


1.44E-03 


2, 


.11 


6.55E-04 


1 


.68 


7.18E-04 


4.93 




81 


3.21E-04 


2, 


.20 


8.24E-04 


-0 


.34 


1.39E-04 


2.41 


mex 


21 


8.67E-03 






2.06E-03 






8.03E-04 






41 


1.87E-03 


2, 


.29 


2.34E-04 


3 


.25 


3.72E-04 


1.15 




81 


3.85E-04 


2, 


.32 


4.93E-05 


2 


.29 


3.82E-05 


3.34 


mic 


21 


5.53E-03 






1.45E-04 






8.15E-04 






41 


1.77E-03 


1, 


.70 


7.00E-05 


1 


.09 


9.21E-04 


-0.18 




81 


3.73E-04 


2, 


.29 


3.45E-04 


-2 


.34 


4.61E-05 


4.40 


plo 


21 


1.04E-03 






1.52E-03 






4.49E-04 






41 


5.06E-04 


1 


.08 


6.05E-04 


1 


.38 


l.OOE-04 


2.24 




81 


7.97E-05 


2 


.71 


1.77E-04 


1 


.81 


2.80E-05 


1.87 


swa 


21 


2.57E-03 






1.08E-03 






1.77E-03 






41 


6.01E-04 


2, 


.17 


3.73E-04 


1 


.58 


4.54E-04 


2.03 




81 


1.50E-04 


2 


.04 


5.45E-05 


2 


.82 


9.28E-05 


2.33 


uch 


21 


1.14E-02 






3.14E-03 






3.05E-03 






41 


1.91E-03 


2, 


.67 


1.08E-04 


5, 


.04 


7.05E-04 


2.19 




81 


3.78E-04 


2, 


.38 


7.68E-04 


-2 


.88 


1.84E-04 


1.97 
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Table 3: Quadratic error for problem 2 



Region 


oize 


S tVllCtllTCd 






UlSblVi CSfla 




UlSblVl CSfl}) 




dom 


21 


9.68E-04 






4.21E-04 




5.40E-04 






41 


2.26E-04 


2, 


.18 


7.54E-05 


2.57 


1.44E-04 


1.98 




81 


4.75E-05 


2 


.29 


1.04E-05 


2.91 


3.16E-05 


2.23 


eng 


21 


9.72E-04 






6.31E-04 




7.68E-04 






41 


1.83E-04 


2, 


.50 


2.49E-05 


4.83 


2.63E-04 


1.60 




81 


5.00E-05 


1 


.90 


4.40E-05 


-0.84 


1.30E-04 


1.04 


hab 


21 


1.08E-03 






1.95E-04 




3.54E-04 






41 


3.90E-04 


1, 


.52 


2.79E-05 


2.91 


6.58E-05 


2.51 




81 


9.72E-05 


2, 


.04 


1.37E-05 


1.04 


1.99E-05 


1.76 


mlQ 


21 


1.30E-03 






2.93E-04 




2.34E-03 






41 


2.75E-04 


2, 


.32 


2.52E-05 


3.66 


1.26E-04 


4.37 




81 


8.46E-05 


1 


.73 


1.33E-04 


-2.45 


2.91E-05 


2.15 


mex 


21 


1.96E-03 






5.39E-04 




3.24E-04 






41 


3.00E-04 


2, 


.80 


8.17E-05 


2.82 


9.53E-05 


1.83 




81 


7.54E-05 


2 


.03 


1.46E-05 


2.53 


4.65E-06 


4.44 


mic 


21 


8.01E-04 






2.78E-05 




1.08E-04 






41 


2.03E-04 


2, 


.05 


5.99E-05 


-1.15 


7.51E-05 


0.55 




81 


6.72E-05 


1 


.62 


3.90E-05 


0.63 


4.65E-06 


4.09 


plo 


21 


2.36E-04 






2.87E-04 




1.13E-04 






41 


9.72E-05 


1 


.33 


1.18E-04 


1.33 


2.26E-05 


2.41 




81 


2.28E-05 


2 


.13 


2.31E-05 


2.40 


4.98E-06 


2.22 


swa 


21 


4.63E-04 






2.35E-04 




3.66E-04 






41 


1.14E-04 


2, 


.10 


7.44E-05 


1.72 


8.81E-05 


2.13 




81 


3.37E-05 


1 


.79 


1.15E-05 


2.75 


2.04E-05 


2.15 


uch 


21 


1.33E-03 






9.61E-04 




4.75E-04 






41 


2.59E-04 


2, 


.44 


2.29E-05 


5.59 


1.30E-04 


1.94 




81 


6.14E-05 


2, 


.12 


1.33E-04 


-2.58 


3.00E-05 


2.15 
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Table 4: Quadratic error for problem 3 



Region 


oize 


*S' tVllCtllTCd 






UlSblVi CSfla 




UlSblVi CSfl}) 




dom 


21 


9.70E-04 






4.22E-04 




5.49E-04 






41 


2.26E-04 


2, 


.18 


7.55E-05 


2.57 


1.47E-04 


1.96 




81 


4.74E-05 


2 


.29 


1.04E-05 


2.91 


3.15E-05 


2.27 


eng 


21 


9.69E-04 






6.33E-04 




6.94E-04 






41 


1.81E-04 


2, 


.50 


2.47E-05 


4.85 


2.60E-04 


1.47 




81 


4.98E-05 


1 


.90 


4.43E-05 


-0.85 


1.37E-04 


0.94 


hab 


21 


1.05E-03 






1.96E-04 




3.55E-04 






41 


3.84E-04 


1, 


.51 


2.76E-05 


2.93 


6.57E-05 


2.52 




81 


9.54E-05 


2, 


.05 


1.37E-05 


1.03 


2.00E-05 


1.75 


mlQ 


21 


1.29E-03 






2.91E-04 




2.26E-03 






41 


2.72E-04 


2, 


.33 


2.46E-05 


3.70 


1.24E-04 


4.34 




81 


8.38E-05 


1 


.73 


1.36E-04 


-2.52 


2.87E-05 


2.15 


mex 


21 


2.01E-03 






5.55E-04 




3.18E-04 






41 


3.04E-04 


2, 


.82 


8.10E-05 


2.88 


9.42E-05 


1.82 




81 


7.32E-05 


2 


.09 


1.45E-05 


2.53 


4.65E-06 


4.42 


mic 


21 


8.13E-04 






2.81E-05 




l.lOE-04 






41 


2.03E-04 


2, 


.07 


7.72E-05 


-1.51 


8.61E-05 


0.37 




81 


6.77E-05 


1 


.61 


4.39E-05 


0.83 


4.96E-06 


4.19 


plo 


21 


2.36E-04 






3.05E-04 




1.14E-04 






41 


9.84E-05 


1 


.31 


1.23E-04 


1.36 


2.28E-05 


2.41 




81 


2.31E-05 


2 


.13 


2.30E-05 


2.46 


5.04E-06 


2.22 


swa 


21 


4.44E-04 






2.44E-04 




3.93E-04 






41 


1.07E-04 


2, 


.12 


7.82E-05 


1.70 


9.46E-05 


2.13 




81 


3.32E-05 


1 


.72 


1.21E-05 


2.74 


2.18E-05 


2.15 


uch 


21 


1.32E-03 






l.OlE-03 




4.97E-04 






41 


2.57E-04 


2, 


.45 


2.29E-05 


5.66 


1.35E-04 


1.95 




81 


6.26E-05 


2 


.07 


1.35E-04 


-2.60 


3.09E-05 


2.17 
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where Ei is the quadratic error associated to the numerical solution calculated 
with a grid with Ui points per side. 

The quadratic error for the grids with 81 points per side for the three problems 
is sketched in figures SI [S] and HI 

Some conclusions can be drawn from the numerical results. The strong non 
convexities on some boundaries are clearly reflected in the error magnitudes 
and in a slight decrease of the empirical order; this can be explained in terms 
of the presence of elongated elements. For instance, in the DistMesha grids, 
the boundary triangles must preserve the boundary shape, so the boundary 
nodes are kept fixed (see fig. (2]) and in consequence some elongated elements 
are produced. Nevertheless, one must also note that the very elongated cells in 
the structured grids have less negative effect in the solution than expected. 
It can also be seen that in a number of problems the solution calculated with 
finite elements is slightly more accurate than than solution obtained with finite 
differences. But the former often required a larger number of unknowns, which 
increases the computational effort required for the calculations. 
However, at the end, one conclusion arises: in the experiments, no method 
seems to be clearly superior. The numerical solutions obtained by second order 
finite differences in the structured grids generated by variational methods are 
essentially as accurate as that obtained by second order finite elements using 
Delaunay-like triangulations. Having in mind this fact, an additional important 
issue must be discussed: the simplicity of use of the data structures required for 
finite differences. Since structured grids are logically rectangular, equation 1^ is 
algorithmically as simple as a nine point discrete laplacian. This is an advantage 
over the more complicated data structures required for triangulations, since the 
numerical results do not show a significative improvement. 



5 Conclusions and future work 

In this paper, we selected irregular planar regions instead of rectangular ones, 
which are often a poor approximation to real-world domains. In addition, the 
side selection was quite arbitrary, and no boundary was processed in any way 
in order to avoid simplifying the problems. Thus, strong non convex bound- 
aries are clearly refiected in the quadratic errors [e.g. Great Britain, M19 and 
Ucha), but these choices were done so because we wanted to deal with a "hard" 
problem, although it must be acknowledged that boundary preprocessing is an 
excellent strategy to improve numerical results in general grids. 
However, the use of irregular fixed boundaries is precisely what leads to the 
main conclusion: as follows from the numerical experiments, no method seems 
to be notably more accurate than the others; difference schemes applied on 
structured meshes can indeed be used to produce reliable approximations to the 
solution of the test problems. In other words, triangulations are not the only 
reliable choice for such regions, it is also possible to generate accurate numerical 
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Figure 4: quadratic error for problem 1 



solutions using grids generated by variational methods in very irregular regions, 
and this approach has not been thoroughly studied yet. As we mentioned, this 
is an important issue, since differences are based on logically rectangular grids, 
and, in consequence, their algorithmic implementations are rather simple. 
Our current research deals with time-dependent partial differential equations 
on irregular domains, and the corresponding results will be reported in a future 
paper. 
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